% function calculate e-predicted coherency
% to be used with multiple coherency calculations
% Caution: alpha only - not tested !!! 29/7
% TESTED WITH GEOTOOLS CORRECT 12pm, 30.07.2002
% e-predicted hxhy-ex
% Reference : Ellinghaus, A. 1997. ProcMT User's manual. Metronix GmbH, Braunschweig, Germany
% latest date 30.07.2002
% Trying to adapt this to compute coherence between Ex(local) and
% Ex(remote)

function[coh] = cohe_h(a,b,data,flag),


% a & b are corresponding impedance (Zxx Zxy) or (Zyx Zyy)
% data is array of cross spectra - SPMmatrix
% 
%flag should be either 'Ex' or 'Ey';

j = sqrt(-1);

HxHx = data(1,1);
HyHy = data(2,2);
HzHz = data(3,3);
ExEx = data(4,4);
EyEy = data(5,5);

HxHy = data(2,1) - j*data(1,2);
HyHx = conj(HxHy);
HxEx = data(4,1) - j*data(1,4);
ExHx = conj(HxEx);
HxEy = data(5,1) - j*data(1,5);
EyHx = conj(HxEy);
ExHy = data(4,2) + j*data(2,4);
HyEx = conj(ExHy);
HyEy = data(5,2) - j*data(2,5);
EyHy = conj(HyEy);


% Ex_p = Zxx*Hx + Zxy*Hy
if flag == 'Ex',
ZpZ = a*HxEx + b*HyEx;
elseif flag == 'Ey',
ZpZ = a*HxEy + b*HyEy;
end;


ZpX = a*HxHx + b*HyHx;
ZpY = a*HxHy + b*HyHy;

%ZpZp = a XZp + b YZp = a (ZpX) + b (ZpY) ;

ZpZp = a*conj(ZpX) + b*conj(ZpY);


if flag == 'Ex',
	ZpZp = ZpZp*ExEx;
elseif flag == 'Ey',
	ZpZp = ZpZp*EyEy;
end;

if abs(ZpZp) > 0,
	Ccoh = ZpZ/sqrt(ZpZp);
	else,
	Ccoh = 1 + j;
end;

coh = abs(Ccoh);

if coh > 1.0,
	coh = 1/coh;
end;

